clear all
version 15
cd "_______"
set more off

local idval=7 // idval=7 for Beijing , idval=67 for Shanghai
forv yearval=2001/2010{
use "./data/pm10data.dta", clear
qui keep if id==`idval' // city by city
qui keep if year==`yearval' // year by year
scalar cff1=0.15
sum pmconc if pmconc<=float(cff1)
gen Fxnp=r(N)/_N  
gen cens=0

*****************
* CENSORED MLE *
*****************
scalar lftcff1=0.135   
scalar rtcff1=0.18
replace cens=1 if pmconc>float(lftcff1) & pmconc <=float(cff1)  
gb2lfit_mod pmconc, cdf(pmparacdf2)  censvar(cens)
local ahat=e(ba)
local bhat=e(bb)
local phat=e(bp)
local qhat=e(bq)
local atx=cff1
gen Fxp_cff=ibeta(`phat',`qhat', (`atx'/`bhat')^`ahat'/(1+(`atx'/`bhat')^`ahat'))
local atx=rtcff1
gen Fxp_rtcff=ibeta(`phat',`qhat', (`atx'/`bhat')^`ahat'/(1+(`atx'/`bhat')^`ahat'))
gen prop=Fxnp-Fxp_cff
gen inrange=Fxp_rtcff-Fxp_cff
gen prop_inrange=prop/inrange

local cityname=city[1] 

cumul pmconc, gen(pmecdf) 
gen cutoff=pmconc*0+cff1
line pmecdf pmconc if pmconc<=0.4, sort xlabel(0(0.1)0.4) lcolor(black) graphregion(color(white)) ///
|| line pmparacdf pmconc if pmconc<=0.4, lcolor(navy) sort ///
|| line pmecdf cutoff, lcolor(black) lpattern(dot) title("`cityname': `yearval'")  ///
name(`cityname'`yearval'_MLE0, replace) aspect(0.8) ///
legend(order(1 "Reported Data" 2 "Estimated Counterfactual Distribution") size(4) region(c(none)) bm(zero)) ytitle("") 

keep prop inrange prop_inrange year 
keep in 1
if `yearval'==2001{
save "./data/`cityname'_inrangeplotestimates.dta"
}
else {
append using "./data/`cityname'_inrangeplotestimates.dta"
save "./data/`cityname'_inrangeplotestimates.dta", replace
}
}


graph combine  Beijing2007_MLE0  Beijing2008_MLE0, ///
rows(1) graphregion(color(white)) ysize(1.5) xsize(3)  
graph export "./outputgraphs/main/Figure4A.pdf", replace


use "./data/Beijing_inrangeplotestimates.dta", replace
graph bar inrange prop prop_inrange, over(year) bar(1, color(gs1)) bar(2, color(gs5)) bar(3, color(gs9)) ///
graphregion(color(white)) title("Beijing: 2001-2010") ///
legend(order(1 "Prop. of Manipulable Data" 2 "Prop. of Manipulated Data" 3 "Prop. of In-range Manipulation") ///
rows(2) size(3) region(c(none)) bm(zero)) ytitle("") 
graph export "./outputgraphs/main/Figure5.pdf", replace
